function gS = quadraticGradient(A, X)

% d (X'*A*X) / d X

[N D] = size(X);

gS = zeros(D, D, N, D);
for i = 1:D
    for j = 1:D
        gS(i, j, :, i) = gS(i, j, :, i) + reshape(A*X(:, j), [1 1 N]);
        gS(i, j, :, j) = gS(i, j, :, j) + reshape(X(:, i)'*A, [1 1 N]);
    end
end
